home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 3: Developer Tools / Linux Cubed Series 3 - Developer Tools.iso / devel / lang / fortran / pgplot5.1 / pgplot5 / pgplot5.1.0 / examples-src / pgdemo3.f < prev    next >
Encoding:
Text File  |  1995-04-17  |  18.9 KB  |  651 lines
  1.       PROGRAM PGDEM3
  2. C-----------------------------------------------------------------------
  3. C Demonstration program for PGPLOT contouring routines.
  4. C-----------------------------------------------------------------------
  5.       INTEGER PGBEG
  6.       WRITE (*,'(A)') ' Demonstration of PGPLOT contouring routines'
  7. C
  8. C Call PGBEG to initiate PGPLOT and open the output device; PGBEG
  9. C will prompt the user to supply the device name and type.
  10. C
  11.       IF (PGBEG(0,'?',1,1) .NE. 1) STOP
  12. C
  13. C Call the demonstration subroutines.
  14. C
  15.       WRITE (*,'(A)') ' Routine PGCONT'
  16.       CALL PGEX31
  17.       WRITE (*,'(A)') ' Routine PGCONS'
  18.       CALL PGEX32
  19.       WRITE (*,'(A)') ' Routine PGCONT with PGCONL labels'
  20.       CALL PGEX36
  21.       WRITE (*,'(A)') ' Routine PGCONB'
  22.       CALL PGEX33
  23.       WRITE (*,'(A)') ' Routine PGCONX with arrow labels'
  24.       CALL PGEX37
  25.       WRITE (*,'(A)') ' Routine PGCONX'
  26.       CALL PGEX34
  27.       WRITE (*,'(A)') ' Routine PGVECT'
  28.       CALL PGEX35
  29. C
  30. C Finally, call PGEND to terminate things properly.
  31. C
  32.       CALL PGEND
  33. C-----------------------------------------------------------------------
  34.       END
  35.  
  36.       SUBROUTINE PGEX31
  37. C-----------------------------------------------------------------------
  38. C Demonstration of contouring routine PGCONT.
  39. C-----------------------------------------------------------------------
  40.       INTEGER I,J
  41.       REAL F(40,40),FMIN,FMAX,ALEV,TR(6)
  42.       DATA TR/0.,1.,0.,0.,0.,1./
  43. C
  44. C Compute a suitable function.
  45. C
  46.       FMIN = 0.0
  47.       FMAX = 0.0
  48.       DO 20 I=1,40
  49.           DO 10 J=1,40
  50.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  51.      1                     (I-J)/40.0
  52.               FMIN = MIN(F(I,J),FMIN)
  53.               FMAX = MAX(F(I,J),FMAX)
  54.    10     CONTINUE
  55.    20 CONTINUE
  56. C
  57. C Clear the screen. Set up window and viewport.
  58. C
  59.       CALL PGPAGE
  60.       CALL PGSVP(0.05,0.95,0.05,0.95)
  61.       CALL PGSWIN(1.0,40.0,1.0,40.0)
  62.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  63.       CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONT')
  64. C
  65. C Draw the map.  PGCONT is called once for each contour, using
  66. C different line attributes to distinguish contour levels.
  67. C
  68.       CALL PGBBUF
  69.       DO 30 I=1,21
  70.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  71.           IF (MOD(I,5).EQ.0) THEN
  72.               CALL PGSLW(3)
  73.           ELSE
  74.               CALL PGSLW(1)
  75.           END IF
  76.           IF (I.LT.10) THEN
  77.               CALL PGSCI(2)
  78.               CALL PGSLS(2)
  79.           ELSE
  80.               CALL PGSCI(3)
  81.               CALL PGSLS(1)
  82.           END IF
  83.           CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR)
  84.    30 CONTINUE
  85.       CALL PGSLW(1)
  86.       CALL PGSLS(1)
  87.       CALL PGSCI(1)
  88.       CALL PGEBUF
  89.       END
  90.  
  91.       SUBROUTINE PGEX32
  92. C-----------------------------------------------------------------------
  93. C Demonstration of contouring routine PGCONS.
  94. C-----------------------------------------------------------------------
  95.       INTEGER I,J
  96.       REAL F(40,40),FMIN,FMAX,ALEV,TR(6)
  97.       DATA TR/0.,1.,0.,0.,0.,1./
  98. C
  99. C Compute a suitable function.
  100. C
  101.       FMIN = 0.0
  102.       FMAX = 0.0
  103.       DO 20 I=1,40
  104.           DO 10 J=1,40
  105.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  106.      1                     (I-J)/40.0
  107.               FMIN = MIN(F(I,J),FMIN)
  108.               FMAX = MAX(F(I,J),FMAX)
  109.    10     CONTINUE
  110.    20 CONTINUE
  111. C
  112. C Clear the screen. Set up window and viewport.
  113. C
  114.       CALL PGPAGE
  115.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  116.       CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONS')
  117. C
  118. C Draw the map.  PGCONS is called once for each contour, using
  119. C different line attributes to distinguish contour levels.
  120. C
  121.       CALL PGBBUF
  122.       DO 40 I=1,21
  123.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  124.           IF (MOD(I,5).EQ.0) THEN
  125.               CALL PGSLW(3)
  126.           ELSE
  127.               CALL PGSLW(1)
  128.           END IF
  129.           IF (I.LT.10) THEN
  130.               CALL PGSCI(2)
  131.               CALL PGSLS(2)
  132.           ELSE
  133.               CALL PGSCI(3)
  134.               CALL PGSLS(1)
  135.           END IF
  136.           CALL PGCONS(F,40,40,1,40,1,40,ALEV,-1,TR)
  137.    40 CONTINUE
  138.       CALL PGSLW(1)
  139.       CALL PGSLS(1)
  140.       CALL PGSCI(1)
  141.       CALL PGEBUF
  142.       END
  143.  
  144.       SUBROUTINE PGEX33
  145. C-----------------------------------------------------------------------
  146. C Demonstration of contouring routine PGCONB.
  147. C-----------------------------------------------------------------------
  148.       REAL BLANK
  149.       PARAMETER (BLANK=-1.2E20)
  150.       INTEGER I,J
  151.       REAL F(40,40),FMIN,FMAX,ALEV,TR(6),X,Y,R
  152.       DATA TR/0.,1.,0.,0.,0.,1./
  153. C
  154. C Compute a suitable function.
  155. C
  156.       FMIN = 0.0
  157.       FMAX = 0.0
  158.       DO 20 I=1,40
  159.           DO 10 J=1,40
  160.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  161.      1                     (I-J)/40.0
  162.               FMIN = MIN(F(I,J),FMIN)
  163.               FMAX = MAX(F(I,J),FMAX)
  164.    10     CONTINUE
  165.    20 CONTINUE
  166. C
  167. C "Blank" the data outside an annulus.
  168. C
  169.       DO 60 I=1,40
  170.           DO 50 J=1,40
  171.               R = SQRT((I-20.5)**2 + (J-20.5)**2)
  172.               IF (R.GT.20.0 .OR. R.LT.3.0) F(I,J) = BLANK
  173.    50     CONTINUE
  174.    60 CONTINUE
  175. C
  176.       CALL PGPAGE
  177.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  178.       CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONB')
  179.       CALL PGBBUF
  180.       DO 80 I=1,21
  181.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  182.           IF (MOD(I,5).EQ.0) THEN
  183.               CALL PGSLW(3)
  184.           ELSE
  185.               CALL PGSLW(1)
  186.           END IF
  187.           IF (I.LT.10) THEN
  188.               CALL PGSCI(2)
  189.               CALL PGSLS(2)
  190.           ELSE
  191.               CALL PGSCI(3)
  192.               CALL PGSLS(1)
  193.           END IF
  194.           CALL PGCONB(F,40,40,1,40,1,40,ALEV,-1,TR,BLANK)
  195.    80 CONTINUE
  196.       CALL PGEBUF
  197. C
  198. C Mark the blanked points for easy identification.
  199. C
  200.       CALL PGBBUF
  201.       CALL PGSCI(1)
  202.       DO 100 I=1,40
  203.           DO 90 J=1,40
  204.               IF (F(I,J).EQ.BLANK) THEN
  205.                   X = TR(1) + REAL(I)*TR(2) + REAL(J)*TR(3)
  206.                   Y = TR(4) + REAL(I)*TR(5) + REAL(J)*TR(6)
  207.                   CALL PGPT(1, X, Y, -1)
  208.               END IF
  209.    90     CONTINUE
  210.   100 CONTINUE
  211.       CALL PGEBUF
  212.       END
  213.  
  214.       SUBROUTINE PGEX34
  215. C-----------------------------------------------------------------------
  216. C This program is intended to demonstrate the use of the PGPLOT routine
  217. C PGCONX. As an example, we take data defined on a sphere. We want to
  218. C draw a contour map of the data on an equal-area projection of the
  219. C surface of the sphere; we choose the Hammer-Aitoff equal-area
  220. C projection centered on Declination (latitude) 0, Right Ascension
  221. C (longitude) 0. The data are defined at 2-degree intervals in both
  222. C coordinates. We thus need a data array dimensioned 181 by 91; the
  223. C array index runs from -90 to +90 in declination (91 elements) and
  224. C from -180 to +180 in right ascension (181 elements). The data at -180 
  225. C and +180 must be identical, of course, but they need to be duplicated 
  226. C in the array as these two longitudes appear on opposite sides of the
  227. C map. 
  228. C-----------------------------------------------------------------------
  229.       REAL  RPDEG
  230.       PARAMETER (RPDEG=3.1415926/180.0)
  231.       INTEGER I, J
  232.       REAL RA, DEC, B, L, XC(181), YC(181)
  233.       REAL Q(181,91), C(9)
  234.       EXTERNAL PLCAIT
  235. C
  236. C Call PGENV to create a rectangular window of 4 x 2 units. This is 
  237. C the bounding rectangle of the plot. The JUST argument is 1
  238. C to get equal scales in x and y.  
  239. C
  240.       CALL PGBBUF
  241.       CALL PGENV(-2.0, 2.0, -1.0, 1.0, 1, -2)
  242.       CALL PGLAB('Right Ascension', 'Declination', ' ')
  243.       CALL PGMTXT('t',2.0,0.0,0.0,
  244.      1           'Contouring on a non-Cartesian grid using PGCONX')
  245.       CALL PGSCH(0.6)
  246.       CALL PGMTXT('b',8.0,0.0,0.0,
  247.      1            'Hammer-Aitoff Equal-Area Projection of the Sphere')
  248.       CALL PGSCH(1.0)
  249. C
  250. C Draw 7 lines of constant longitude at longitude 0, 60, 120, ..., 
  251. C 360 degrees. Each line is made up of 90 straight-line segments.
  252. C
  253.       DO 20 J=1,7
  254.           RA = (-180.+(J-1)*60.)*RPDEG
  255.           DO 10 I=1,91
  256.               DEC = 2*(I-46)*RPDEG
  257.               CALL AITOFF(DEC,RA,XC(I),YC(I))
  258.    10     CONTINUE
  259.           CALL PGLINE(91,XC,YC)
  260.    20 CONTINUE
  261. C
  262. C Draw 5 lines of constant latitude at latitudes -60, -30, 0, 30, 
  263. C 60 degrees. Each line is made up of 360 straight-line segments.
  264. C
  265.       DO 40 J=1,5
  266.           DEC = (-60.+(J-1)*30.)*RPDEG
  267.           DO 30 I=1,181
  268.               RA = 2*(I-91)*RPDEG
  269.               CALL AITOFF(DEC,RA,XC(I),YC(I))
  270.    30     CONTINUE
  271.           CALL PGLINE(181,XC,YC)
  272.    40 CONTINUE
  273.       CALL PGEBUF
  274. C
  275. C Compute the data to be contoured. In practice the data might be read
  276. C in from an external file. In this example the data are computed: they
  277. C are the galactic latitudes of the points on the sphere. Thus the 
  278. C contours will be lines of constant galactic latitude.
  279. C
  280.       DO 60 J=1,91
  281.           DEC = 2*(J-46)*RPDEG
  282.           DO 50 I=1,181
  283.               RA = 2*(I-91)*RPDEG
  284.               CALL GALACT(RA, DEC, B,L)
  285.               Q(I,J) = B
  286.    50     CONTINUE
  287.    60 CONTINUE
  288. C
  289. C Draw the contour map using PGCONX. Contours at 0, 20, 40, 60, 80.
  290. C
  291.       DO 70 I=1,9
  292.           C(I) = -100.0 +I*20.0
  293.    70 CONTINUE
  294.       CALL PGBBUF
  295.       CALL PGSCI(2)
  296.       CALL PGCONX(Q, 181, 91, 1, 181, 1, 91, C, 9, PLCAIT)
  297.       CALL PGSCI(1)
  298.       CALL PGEBUF
  299.       END
  300.  
  301.       SUBROUTINE PLCAIT(VISBLE, X, Y, Z)
  302.       INTEGER VISBLE
  303.       REAL X,Y,Z
  304. C-----------------------------------------------------------------------
  305. C Plotting subroutine for PGCONX. This routine converts the input
  306. C coordinates (latitude and longitude) into the projected coordinates
  307. C (x and y), and moves or draws as requested by VISBLE.
  308. C-----------------------------------------------------------------------
  309.       REAL  RPDEG
  310.       PARAMETER (RPDEG=3.1415926/180.0)
  311.       REAL B, L, XWORLD, YWORLD
  312.       B = 2.0*(Y-46.0)*RPDEG
  313.       L = 2.0*(X-91.0)*RPDEG
  314.       CALL AITOFF(B, L, XWORLD, YWORLD)
  315.       IF (VISBLE.EQ.0) THEN
  316.           CALL PGMOVE(XWORLD, YWORLD)
  317.       ELSE
  318.           CALL PGDRAW(XWORLD, YWORLD)
  319.       END IF
  320.       END
  321.  
  322.       SUBROUTINE AITOFF(B,L,X,Y)
  323. C-----------------------------------------------------------------------
  324. C Hammer-Aitoff projection.
  325. C
  326. C       Input: latitude and longitude (B,L) in radians
  327. C       Output: cartesian (X,Y) in range +/-2, +/-1
  328. C-----------------------------------------------------------------------
  329.       REAL L,B,X,Y,L2,DEN
  330. C
  331.       L2 = L/2.0
  332.       DEN = SQRT(1.0+COS(B)*COS(L2))
  333.       X = 2.0*COS(B)*SIN(L2)/DEN
  334.       Y = SIN(B)/DEN
  335.       END
  336.  
  337.       SUBROUTINE GALACT(RA,DEC,GLAT,GLONG)
  338. C-----------------------------------------------------------------------
  339. C Convert 1950.0 equatorial coordinates (RA, DEC) to galactic
  340. C latitude and longitude (GLAT, GLONG).
  341. C
  342. C Arguments:
  343. C  RA, DEC (input): 1950.0 RA and Dec (radians).
  344. C  GLAT, GLONG (output): galactic latitude and longitude 
  345. C      (degrees).
  346. C
  347. C Reference: e.g., D. R. H. Johnson and D. R. Soderblom, A. J. v93, 
  348. C  p864 (1987).
  349. C-----------------------------------------------------------------------
  350.       REAL RA, RRA, DEC, RDEC, CDEC, R(3,3), E(3), G(3)
  351.       REAL RADDEG, GLAT, GLONG
  352.       INTEGER I, J
  353.       DATA R/-.066988740D0, .492728466D0,-.867600811D0,-.872755766D0,
  354.      $       -.450346958D0,-.188374601D0,-.483538915D0, .744584633D0,
  355.      $        .460199785D0/
  356.       DATA RADDEG/57.29577951D0/
  357. C-----------------------------------------------------------------------
  358.       RRA = RA
  359.       RDEC = DEC
  360.       CDEC = COS(RDEC)
  361.       E(1) = CDEC*COS(RRA)
  362.       E(2) = CDEC*SIN(RRA)
  363.       E(3) = SIN(RDEC)
  364.       DO 20 I=1,3
  365.           G(I) = 0.0
  366.           DO 10 J=1,3
  367.               G(I) = G(I) + E(J)*R(I,J)
  368.    10     CONTINUE
  369.    20 CONTINUE
  370.       GLAT  = ASIN(G(3))*RADDEG
  371.       GLONG = ATAN2(G(2),G(1))*RADDEG
  372.       IF (GLONG.LT.0.0) GLONG = GLONG+360.0
  373.       RETURN
  374. C-----------------------------------------------------------------------
  375.       END
  376.  
  377.       SUBROUTINE PGEX35
  378. C-----------------------------------------------------------------------
  379. C Program to demonstrate the use of PGVECT along with
  380. C PGCONB by illustrating the flow around a cylinder with circulation.
  381. C
  382. C          NX      total # of axial stations
  383. C          NY      total # of grid pts in y (or r) direction
  384. C-----------------------------------------------------------------------
  385.       INTEGER MAXX, MAXY
  386.       PARAMETER (MAXX=101,MAXY=201)
  387.       INTEGER NX, NY, I, J, NC
  388.       REAL PI, A, GAMMA, VINF, XMAX, XMIN, YMAX, YMIN, DX, DY
  389.       REAL CPMIN, R2, BLANK
  390.       REAL CP(MAXX,MAXY),X(MAXX),Y(MAXY),U(MAXX,MAXY),V(MAXX,MAXY),
  391.      1   PSI(MAXX,MAXY)
  392.       REAL TR(6),C(10)
  393.       DATA BLANK/-1.E10/
  394. C
  395. C compute the flow about a cylinder with circulation
  396. C
  397. C define various quantities
  398. C
  399. C pi
  400.       PI = ACOS(-1.)
  401. C number of points in the x and y directions
  402.       NX = 31
  403.       NY = 31
  404. C cylinder radius
  405.       A = 1.
  406. C circulation strength
  407.       GAMMA = 2.
  408. C freestream velocity
  409.       VINF = 1.
  410. C max and min x and y
  411.       XMAX = 3.*A
  412.       XMIN = -3.*A
  413.       YMAX = 3.*A
  414.       YMIN = -3.*A
  415. C point spacing
  416.       DX = (XMAX-XMIN)/(NX-1)
  417.       DY = (YMAX-YMIN)/(NY-1)
  418. C compute the stream function, Cp, and u and v velocities
  419.       CPMIN =1.E10
  420.       DO 20 I=1,NX
  421.          X(I) = XMIN+DX*(I-1)
  422.          DO 10 J=1,NY
  423.             Y(J) = YMIN+DY*(J-1)
  424.             R2 = X(I)**2+Y(J)**2
  425.             IF (R2.GT.0.) THEN
  426.                PSI(I,J) = VINF*Y(J)*(1.-A**2/R2)
  427.      1            +GAMMA/(2.*PI)*0.5*ALOG(R2/A)
  428.                U(I,J) = VINF*(1.+A**2/R2-2.*A**2*X(I)**2/R2**2)
  429.      1            +GAMMA/(2.*PI)*Y(J)/R2
  430.                V(I,J) = VINF*X(I)*(-2.*A**2*Y(J)/R2**2)
  431.      1            +GAMMA/(2.*PI)*X(I)/R2
  432.                CP(I,J) = 1.-(U(I,J)**2+V(I,J)**2)/VINF**2
  433.             ELSE
  434.                PSI(I,J) = 0.
  435.                U(I,J) = 0.
  436.                V(I,J) = 0.
  437.                CP(I,J) = 0.
  438.             END IF
  439.             IF (R2.LT.A**2) THEN
  440.                U(I,J) = BLANK
  441.                V(I,J) = BLANK
  442.             ELSE
  443.                CPMIN = MIN(CPMIN,CP(I,J))
  444.             END IF
  445.    10    CONTINUE
  446.    20 CONTINUE
  447. C
  448. C grid to world transformation
  449. C
  450.       TR(1)=X(1)-DX
  451.       TR(2)=DX
  452.       TR(3)=0.0
  453.       TR(4)=Y(1)-DY
  454.       TR(5)=0.0
  455.       TR(6)=DY
  456. C
  457.       CALL PGENV (X(1),X(NX),Y(1),Y(NY),1,0)
  458.       CALL PGIDEN
  459.       CALL PGLAB ('X','Y','Flow About a Cylinder with Circulation')
  460. C
  461. C contour plot of the stream function (streamlines)
  462. C
  463.       NC=5
  464.       C(1)=1.
  465.       C(2)=.5
  466.       C(3)=0.
  467.       C(4)=-.5
  468.       C(5)=-1.
  469.       CALL PGCONT (PSI,MAXX,MAXY,1,NX,1,NY,C,NC,TR)
  470. C
  471. C draw cylinder
  472. C
  473.       CALL PGBBUF
  474.       CALL PGSCI (0)
  475.       CALL PGSFS (1)
  476.       CALL PGCIRC (0.,0.,A*1.1)
  477.       CALL PGSFS (2)
  478.       CALL PGSCI (14)
  479.       CALL PGCIRC (0.0, 0., A)
  480.       CALL PGSCI (1)
  481.       CALL PGEBUF
  482. C
  483. C vector plot
  484. C
  485.       CALL PGSAH (2, 45.0, 0.7)
  486.       CALL PGSCH (0.3)
  487.       CALL PGVECT (U,V,MAXX,MAXY,2,NX-1,2,NY-1,0.0,0,TR,-1.E10)
  488.       CALL PGSCH(1.0)
  489. C
  490. C finished
  491. C
  492.       RETURN
  493. C----------------------------------------------------------------------
  494.       END
  495.  
  496.       SUBROUTINE PGEX37
  497. C-----------------------------------------------------------------------
  498. C Demonstration of contouring routine PGCONX.
  499. C-----------------------------------------------------------------------
  500.       INTEGER I,J
  501.       REAL F(40,40),FMIN,FMAX,ALEV
  502.       EXTERNAL PLCARO
  503. C
  504. C Compute a suitable function.
  505. C
  506.       FMIN = 0.0
  507.       FMAX = 0.0
  508.       DO 20 I=1,40
  509.           DO 10 J=1,40
  510.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  511.      1                     (I-J)/40.0
  512.               FMIN = MIN(F(I,J),FMIN)
  513.               FMAX = MAX(F(I,J),FMAX)
  514.    10     CONTINUE
  515.    20 CONTINUE
  516. C
  517. C Clear the screen. Set up window and viewport.
  518. C
  519.       CALL PGPAGE
  520.       CALL PGSVP(0.05,0.95,0.05,0.95)
  521.       CALL PGSWIN(1.0,40.0,1.0,40.0)
  522.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  523.       CALL PGMTXT('t',1.0,0.0,0.0,
  524.      :            'Contouring using PGCONX with arrows')
  525. C
  526. C Draw the map.  PGCONX is called once for each contour, using
  527. C different line attributes to distinguish contour levels.
  528. C
  529.       CALL PGBBUF
  530.       DO 30 I=1,21
  531.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  532.           IF (MOD(I,5).EQ.0) THEN
  533.               CALL PGSLW(3)
  534.           ELSE
  535.               CALL PGSLW(1)
  536.           END IF
  537.           IF (I.LT.10) THEN
  538.               CALL PGSCI(2)
  539.               CALL PGSLS(2)
  540.           ELSE
  541.               CALL PGSCI(3)
  542.               CALL PGSLS(1)
  543.           END IF
  544.           CALL PGCONX(F,40,40,1,40,1,40,ALEV,-1,PLCARO)
  545.    30 CONTINUE
  546.       CALL PGSLW(1)
  547.       CALL PGSLS(1)
  548.       CALL PGSCI(1)
  549.       CALL PGEBUF
  550.       END
  551.  
  552.       SUBROUTINE PGEX36
  553. C-----------------------------------------------------------------------
  554. C Demonstration of contouring routine PGCONT and PGCONL.
  555. C-----------------------------------------------------------------------
  556.       INTEGER I,J
  557.       REAL F(40,40),FMIN,FMAX,ALEV,TR(6)
  558.       CHARACTER*32 LABEL
  559.       DATA TR /0.0, 1.0, 0.0, 0.0, 0.0, 1.0/
  560. C
  561. C Compute a suitable function.
  562. C
  563.       FMIN = 0.0
  564.       FMAX = 0.0
  565.       DO 20 I=1,40
  566.           DO 10 J=1,40
  567.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  568.      1                     (I-J)/40.0
  569.               FMIN = MIN(F(I,J),FMIN)
  570.               FMAX = MAX(F(I,J),FMAX)
  571.    10     CONTINUE
  572.    20 CONTINUE
  573. C
  574. C Clear the screen. Set up window and viewport.
  575. C
  576.       CALL PGPAGE
  577.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  578.       CALL PGMTXT('t',1.0,0.0,0.0,
  579.      1            'Contouring using PGCONT and PGCONL labels')
  580. C
  581. C Draw the map.  PGCONT is called once for each contour, using
  582. C different line attributes to distinguish contour levels.
  583. C
  584.       CALL PGBBUF
  585.       DO 40 I=1,21
  586.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  587.           IF (MOD(I,5).EQ.0) THEN
  588.               CALL PGSLW(3)
  589.           ELSE
  590.               CALL PGSLW(1)
  591.           END IF
  592.           IF (I.LT.10) THEN
  593.               CALL PGSCI(2)
  594.               CALL PGSLS(2)
  595.           ELSE
  596.               CALL PGSCI(3)
  597.               CALL PGSLS(1)
  598.           END IF
  599.           CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR)
  600.    40 CONTINUE
  601.       CALL PGSLW(1)
  602.       CALL PGSLS(1)
  603.       CALL PGEBUF
  604. C
  605. C Label the contours with PGCONL. Only even-numbered contours
  606. C are labelled.
  607. C
  608.       CALL PGBBUF
  609.       DO 50 I=2,21,2
  610.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  611.           WRITE (LABEL,'(I2)') I
  612. C         WRITE (LABEL,'(F8.2)') ALEV
  613.           IF (I.LT.10) THEN
  614.               CALL PGSCI(2)
  615.           ELSE
  616.               CALL PGSCI(3)
  617.           END IF
  618.           CALL PGCONL(F,40,40,1,40,1,40,ALEV,TR,LABEL,16,8)
  619.  50   CONTINUE
  620.       CALL PGSCI(1)
  621.       CALL PGEBUF
  622.       END
  623.  
  624.       SUBROUTINE PLCARO(VISBLE, X, Y, Z)
  625.       INTEGER VISBLE
  626.       REAL X,Y,Z
  627. C-----------------------------------------------------------------------
  628. C Plotting subroutine for PGCONX. This routine labels contours with
  629. C arrows. Arrows point clockwise around minima, anticlockwise around
  630. C maxima. Arrows are drawn on 1/16 of contour line segments.
  631. C-----------------------------------------------------------------------
  632.       REAL XP, YP
  633.       INTEGER I
  634.       SAVE I
  635.       DATA I /0/
  636. C
  637.       I = MOD(I+1,16)
  638.       IF (VISBLE.EQ.0) THEN
  639.           I = 0
  640.           CALL PGMOVE(X, Y)
  641.       ELSE IF (I.EQ.8) THEN
  642. C         -- Draw line segment with arrow at midpoint
  643.           CALL PGQPOS(XP,YP)
  644.           CALL PGARRO(XP, YP, (X+XP)*0.5, (Y+YP)*0.5)
  645.           CALL PGDRAW(X, Y)
  646.       ELSE
  647. C         -- Draw plain line segment
  648.           CALL PGDRAW(X, Y)
  649.       END IF
  650.       END
  651.